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Abstract: We have developed a general formalism to compute Sommerfeld enhancement (SE) 
factors for a multi-state system of fermions, in all possible spin configurations and with generic 
long-range interactions. We show how to include such SE effects in an accurate calculation of the 
thermal relic density for WIMP dark matter candidates. We apply the method to the MSSM and 
perform a numerical study of the relic abundance of neutralinos with arbitrary composition and 
including the SE due to the exchange of the W and Z bosons, photons and Higgses. We find that 
the relic density can be suppressed by a factor of a few in a sizable region of the parameter space, 
mostly for Wino-like neutralino with mass of a few TeV, and up to an order of magnitude close to 
a resonance. 
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1. Introduction 

In the latest years a number of experiments have contributed to develop a deep understanding of 
the standard model of cosmology. Most cosmological parameters have been measured with a level 
of precision hardly foreseeable a decade ago; among these, there is the cold dark matter (CDM) 
contribution to the Universe energy density. E.g., a recent analysis, within the 6-parameter ACDM 
model, of the 7- year WMAP data , combined with the distance measurements from the baryon 
acoustic oscillations in the distribution of galaxies and the recent redetermination of the Hubble 
Constant with the Hubble Space Telescope gives [Q: 



where indicates the ratio between mean density and critical density, and h is the Hubble constant 
in units of 100 kms~^ Mpc^^. 

The nature of the CDM term is still unknown; one of the most attractive scenarios is that dark 
matter is a thermal relic from the early Universe: in this context, stable weakly interacting massive 
particles (WIMPs) are ideal CDM candidates, as their thermal relic abundance is naturally of the 
order of the measured one (for a recent review on WIMP DM see, e.g., jsj). Since the accuracy 
in the experimental determination of ficDM has reached the few per cent level, it is necessary to 
have equally accurate theoretical predictions for the relic abundance of WIMPs. Indeed, there 
may be a number of effects involved in this calculation, requiring a proper treatment and possibly 
making the final result differ even by orders of magnitude from the simple and naive estimate 



ilcDuh^ = 0.1123 ±0.0035 



(68% CL uncerinities) , 



(1.1) 
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that the rehc density scales with the inverse of the total WIMP pair-annihilation cross section in 
the non-relativistic regime, and it takes the appropriate value when the annihilation rate is about 
3 • 10~^^ cm'^ s~'^. Such features include, e.g., threshold or resonance effects, co-annihilation effects 
(namely the interplay among several particles, nearly degenerate in mass, rather than a single 
WIMP) or three-body final states. Lately, some attention has been dedicated the case of WIMPs 
with annihilation cross sections strongly dependent on the velocity of the incoming pair. Such 
interest has been triggered by the fact that cross sections much larger the standard value seems to 
be needed to provide a DM positron source accounting for the rise in the positron fraction in the 
local cosmic-rays that has been measured by PAMELA |^ : the idea is to explain such mismatch in 
terms of the mismatch between the velocity of the WIMPs at freeze-out in the early Universe, soon 
after entering the non-relativistic regime, i.e. for about v/c ^ 0.3, with the velocity of WIMPs in 
our Galactic halo, i.e. about w/c ~ 10""^, see, e.g., 0. A feature of this kind is predicted by to 
the so-called Sommerfeld enhancement namely the increase in the cross section occurring for 
highly non-relativistic particles in the presence of a long-range attractive force that deforms the 
wave function of the incoming pair with respect to the standard plane wave approximation. 

The Sommerfeld enhancement effect in the context of WIMP DM has been studied first in 
Refs. [|[ |l^; these authors considered, within a supersymmetric extension of the standard model of 
elementary particles, the case of pure Wino or pure Higgsino DM and showed that, when these states 
are very massive, at the TeV scale or heavier, the weak force can play the role of a long-range force 
since its carriers, the W and Z bosons, are much lighter. The impact on the relic density calculation 
for this case was discussed in Ref. JT]| . Later, the Sommerfeld enhancement and its implication on 
the relic density has been studied in the context of the minimal dark matter model still 
referring to standard interactions and very massive WIMPs. More recently, in light of the PAMELA 
anomaly, several authors, including, e.g., § ||, |6[ ||, ||, |l|, ||, H, ||, H, have 
considered a few aspects of the connection with WIMP DM, considering models for which a new 
dark force, carried by some yet-to-be-discovered light boson, provides a very large enhancement 
in the cross section. In our analysis we reconsider the case of supersymmetric DM, and discuss 
the effect in the general context of the minimal supersymmetric extension to the standard model 
(MSSM), with arbitrary Wino and Higgsino components in the lightest neutralino,^ applying and 
developing the formalism introduced in Ref. (the derivation of the Sommerfeld enhancements 



from field theory was discussed also in 1 27 , see also ) . Such formalism is more suitable to treat a 
realistic case of neutralino DM than the approach based on implementing a non-relativistic effective 
action followed by previous work on this subject 1^ (see also ||l2|). 

The new calculation that we present in this paper is then used to refine the computation of 
the thermal relic density of neutralinos. To solve the corresponding set of Boltzmann equations we 
use the very accurate method given in the public available DarkSUSY package |2^, appropriately 
modified to include Sommerfeld enhancement effects. Hence, we provide here a tool for very high 
precision estimates of relic abundances of neutralino DM in the MSSM, at a level comparable or 
better than the accuracy in the value for flcuuh^ from current and upcoming observations. 

The paper is organized as follows: In Section 2 we introduce the Sommerfeld enhancement and 
sketch the limits in which it may be relevant, while Section 3 discusses our approach to the case of 
TV-coupled states. In Section 4 we detail the technique we use to implement the Sommerfeld effect 
in the relic density computation. Section 5 discusses our particle physics setup, while in Section 6 
we present the results within this framework. Section 7 concludes. 



^We will restrict to the case of negligible bino fraction since there is no Sommerfeld enhancement on trivial 
representation of 5f7{2)£. 
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2. Sommerfeld enhancement 



Wc start with an short overview of the Sommerfeld effect itself and recap of the literature on 
the subject. The Sommerfeld enhancement is a non-relativistic effect resulting in correcting the 
cross section (in our case of interest the annihilation cross section) due to presence of some "large 
distance force" between the particles in the incoming state. It is generally described as an effect 
of distorting the initial wave function ijj{r) of the incoming two-particle state by a non-relativistic 
potential. This potential is taken to be Yukawa or Coulomb, as the force arises due to exchange of 
massive or massless boson. The whole effect can then be encoded into the ratio between the wave 
function of the incoming, free particle (at r oo) and the distorted one at the point of annihilation 
(at r = 0). Thus one usually defines the enhancement factor being 

which multiplies the cross section 

CfuU = S ■ ao- (2.2) 

Many authors have discussed this effect introducing a new interaction, mediated by some scalar or 
vector boson <j>. If one assumes one species of annihilating particle x with mass much larger than 
the mass of force carrier, ^ m^, and a coupling of the order of the weak scale or larger, then 
one can get a large enhancement of the cross section. This can have a very important implications 
for indirect dark matter searches and relic density computations. 

The standard approach to estimate the enhancement is to introduce an interaction potential in 
the form: 

V = a , (2.3) 

r 

and solve the Schordinger equation for the incoming two-particle state to find the wave function 
distortion; there are however a few issues one should be careful about. First, the coefficient a in the 
potential depends on the nature of the incoming state and the possible type(s) of interaction. In 

particular, when the effect involves fermions, it is different if the interacting particles are Dirac or 
Majorana and whether they are identical or not; this is especially important when co-annihilations 
enter the computation of the relic density. 

To stress other delicate points, it is useful to recap first what are the conditions for the en- 
hancement to be sizable. The Sommerfeld enhancement can be viewed as occurring due to forming 
a loosely bound state due to long-range interaction. In order to have such a bound state the char- 
acteristic Bohr energy of the interaction need to be larger than the kinetic energy. In the limit 
0, this gives a condition a^m^ > m^v^, i.e.: 



< 



a . (2.4) 



For a typical WIMP the coupling is of order a ~ 0.03, so that there can be some sizable enhancement 
only long after freeze-out (which happens for v ~ 0.3). However, if there exists a slightly heavier 
state, then it may happen that just after freeze-out DM particles have enough energy to produce 
it nearly on-shell. At threshold these heavier states are produced with, roughly speaking, zero 
velocity. As we will see later on, if the mass splitting of the DM and the heavier state is small 
enough, this may give rise to important changes in the relic density. 

A second condition comes from the comparison of the range of the Yukawa potential with the 
Bohr radius. In order for the interaction to distort the wave function significantly, the range of the 
potential cannot be much smaller than the Bohr radius of the two-particle state, 

^>^. (2.5) 

mij, arrix 
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In case of very large enhancements, this condition needs to be even stronger, i.e. the range of 
potential has to be much larger than the Bohr radius. However, even in case of enhancements of 
order unity one should treat carefully the regime when w Q^™x- 

When considering a system of two states with small mass splitting dm interacting off-diagonally 
there is another important constraint. If dm is significantly larger than the kinetic energy, it may 
seem that the heavier state cannot be produced, and hence there is no enhancement. However, if 
the potential is strong enough, there still may be an effect, coming from producing the heavier state 
at small distances, where the potential energy is large. Thus the condition reads: 

2Sm<a^m^+£, (2.6) 

meaning that the characteristic Bohr energy of the potential plus the kinetic energy £ is large 
enough to produce the heavier state. Moreover, when dealing with multi-state systems, the picture 
of the Sommerfeld enhancement as an effect of a static, long range force is no longer applicable. This 
is because the exchange of leads to a momentum and energy transfer due to the mass splitting 
dm and one may need to take into account terms of order 0{5m/m) modifying the interaction 
potential. 

When dealing with a specific particle physics setup like, e.g., the MSSM, several such compli- 
cations may intervene at the same time: a simple parametric description is not possible and for a 
proper estimate of the Sommerfeld effect and its impact on the relic density, a computation within 
a fully general formalism is needed. We introduce this in the next Section. 

3. Sommerfeld effect for N coupled states 

We are interested in the general case with N two-particle fermionic states coupled together. They 
interact with "long-range forces" due to the exchange of some boson </>, denoting generically a 
vector, an axial vector or a scalar boson (in the MSSM those correspond to W^, 7, H12 and 
H^). The Sommerfeld enhancement corresponds to computing in the non-relativistic limit the sum 
of ladder diagrams as presented on the Fig. Ilia. The structure of those diagrams can be in general 



a) 




Figure 1: Ladder diagrams for the Sommerfeld enhancement, a) incoming XiXj particles interact with 
exchange of (/)'s (in general different), which can be a scalar, vector or axial vector bosons. In the ladder 
a virtual states Xi'Xj' can be produced and the final annihilation proceed in the channel which can be 
different than initial one. Filled blob represents full annihilation process with any number of SM particles 
in final state, while the empty one its tree level counterpart, b) the same but written in a recursive form; 
sum is over all possible intermediate states and exchanged bosons. 

very complicated, since all states can appear in one diagram, through different interactions. The 
approach we follow to address the this problem is a generalization of the one developed in . We 
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start from writing the recurrence relation for the annihilation amplitudes as presented pictorially 
on Fig. 0b. Let's consider a process of the type 

XaXb XiXi x'iX'j ^ ■ . . ^ SM final states , (3.1) 

where the intermediate pairs XiXj can be the same or different as the initial pair XaXb- The spin of 
the initial pair, which in the non-relativistic limit is a conserved quantity, can be in general either 
in the singlet {S — 0) or triplet (5 — 1) state. For every possible XiXj P^'ir ''^6 g^t a recurrence 
relation for its annihilation amplitude. Denoting this amplitude by Aij and its tree level value by 
A'^j, one obtains in the non-relativistic limit pof : 

A., ip) = 4 ip) - Y: N.,,,, f f^^^ ^ ^"'^'^ , (3.2) 

■' ^ ^ 2m,. 

where the sum is over different Xi'Xj' iiitermediate states and different interactions. Here£ = /2m'f' 
is kinetic energy of incoming pair (at infinity) , with m"'' the reduced mass and p the CM three- 
momentum, 26mij — nii + rrij — {ma + rrib) the mass splitting, gu'^ and gj'j,), are the coupling 
constants and Nij^i/ji is the term containing normalization and combinatorial factors. 

To rewrite the expression above in the form of the Schrodinger equations we make following 
redefinitions: 

= ( A - ^ + 2'^"^^^) (P)' (3-3) 



U°ir)= J d^pe^P ^'A^ip^Po), (3.4) 
^.,(rO = / dpe'P-%jip), (3.5) 



which allows us to rewrite Eq. (|3.2| ) as a differential equation: 

where (f) refers to the particle being exchanged (scalar, vector or axial vector boson). The potential 
has the form 

with Cij^i'j' {(f>) being coefficients depending on the couplings and states involved. An efficient way 
of computing them is explained in Appendix ^ and the results in case of a system involving one 
spin 1/2 Dirac fermion and/or two different Majorana spin 1/2 fermions are summarized in Tab. |l|. 
Whether the potential is attractive or repulsive is hidden in the sign of the coefficients c and depends 
on the interaction type. The exchange of scalars is always attractive in the spin singlet (i.e. overall 
plus sign), but can also be repulsive in the triplet. Vector and axial bosons can give attractive or 
repulsive forces, depending on the charges. 

Notably in the approach outlined here we can split different interaction types (i.e. mediated 
by different bosons or with different couplings etc.) within a ladder diagram and consider them 
separately. The trade-off is that every possible intermediate two-particle state will lead to one 
equation, so in such a general case we will need to consider a set of coupled Schrodinger equations. 
These equations are inhomogeneous; however, since the Sommerfeld enhancement factorizes out 
from the annihilation matrix, as it does not depend on the final states but enters only as a distortion 
of the incoming wave function, one can compute it by solving the associated homogeneous, using 
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Spin singlet 



cj> : 


scalar (F = IL) 


vector (i =70) 


axial (i = 7i75j 


c+-,+- 




g' 




c++,++ 


g' 


-g' 


-35^ 


C",H — 


V2\g,+ f 


V2I5.+P 


-3^|ff,+ P 


Cij,H — 


2Re(5,+5*+) 


2Re(3i+5j+) 


-6Re(g,+5f*+) 




^\g^j\^ + gl + gtj 




-3(2|ffi,f 












ISi+l^ + 23i»3 




-3|gi+|2 - 6p,,g 






-gr+g'j+ - "igiimigij) 


-3ffi+5;+ -6ffRe (gij) 









-I25S 




4\/2giiRe (gij) 





-12\/25iiRe (sij) 




Spin triplet 










c+-,+- 




g' 




c++,++ 




-g' 


9' 


Cii, + - 













2ilm {g*+gj+) 


2ilm ig*+gj+) 


2ilm {g*+gj+) 















-(2l5ij P + + ff?/) + 4ff«,to 


01 |2 2 *2 


-(2lfef + 5?,- + .gr/) + 45f,,5jj 


C+i,+i 


+2ffffi» 


\g^+? 


+ 255ii 




+2ffRe(ffij) 


g^+g']+ - 2gilm (gij) 


+ 2gRe(gij) 













Cij,ii 











Couplings: 


sSxjrxi^ (+/I.C. ifr j / j), 


fff+T/irx!^ +h.c., 


g^i}Vip(j), where r = l, 7o,7i75 



Table 1: List of all possible coefficients in the potential V^.^,-,{r) for the Sommerfeld effect computation. 
The table includes any annihilation process involving one spin 1/2 Dirac fermion (denoted by + or — 
depending on whether it is a particle or antiparticle) and/or two different Majorana spin f/2 ferrnions 
(denoted by i and j), for an even partial wave. Couplings are defined in the last line for each F, where 
X is a Majorana fermion, tp a Dirac field, (j) is the exchanged boson. When applied to the MSSM, where 
one needs to consider a chargino and one or two neutralinos, some of the coefficients vanish due to the CP 
conservation and couplings of type ga are negligible. The overall " +" sign refers to an attractive force while 
"— " to a repulsive force. Note also that Cij^ki — cli ^j. 



the partial waves decomposition, as described in [ p6| . We will be interested only in the s-wave, but 
it is straightforward (though not always easy numerically) to extend the analysis to higher partial 
waves. 

To reduce these equations in a form more suitable for numerical calculations, after the partial 
wave decomposition we define the reduced radial wave function (f{x) as 

</W = ^P^^' (3-8) 

where N is some normalization constant and p is the value of the CM three-momentum for one of 
the incoming particles (a or b). Since we restrict to the s-wave case I = 0, we will drop the I index 
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from now on. From ( |3.6| ) one get then set of equations^ for the ip^^ 



= 0. (3.9) 



To obtain the enhancement we need to solve this set of equations with appropriate boundary 
conditions. In a: = they are set by the requirement that the solution is regular. In x — > oo 
the solution has to describe one incoming XaXb state and all the possible XiXj states that can be 
produced in the ladder. For the latter there can be two cases: 

1. 25mij < £ - there is enough energy to produce on-shell states XiXj^ 

2. 25mij > £ - there is not enough energy; states XiXj ^'^^ off-shell. 
The radial wave functions behave at infinity as: 

R^\r)^^'- 1^ , (3.10) 

2i r 2i r 

for the incoming pair and 

^^^M-|f "5:.''°""''''' (3.11) 

1 £2^^1211. if off -shell, 

for every other intermediate state XiXj- We use the normalization of the wave function to be such, 
that for the non-interacting case Rab ~ s\n{kab'r)/r (for details see [^). After changing the variable 
to a; = kab • I" and defining qij — rri^J /rn'f' • (1 — 2Smij /£) we get set of boundary conditions for the 
reduced wave functions at x — >■ 00: 

lip"'' - d^ifi"'' = -e-*^ (3.12) 
'i'y/'lijf'''' ~ dx(p^-' =0 if on — shell, 



y^</j''J + da^ip'^ if ofl? - shell. ^^'^^^ 
To check our numerics we can use the unitarity condition, saying that: 

l = |CfP+VVg-|C^^p. (3.14) 



Since the set of boundary conditions depends on the initial particles masses and energy, the solutions 
to the Schrodinger equations, i.e. the wave functions, depend on the initial conditions of the 
incoming pair XaXb' to be precise we should call them then ip^f^^f^^^x) . After solving ( ^.g| ) the 
(co-)annihilation cross section of the pair XaXb is determined, up to the kinematical factor, by 

'^(ab)^T.€b)-K\"' (3-15) 

where the enhancement factors with our normalization are 

S'(Lb)-\d.^U\"^o- (3.16) 

It is important to note, that the computation of the enhancement depends on the spin state of 
the initial 2-body state (see Tab. |l|). This means that one needs to project each annihilation cross 
section into two parts, one for the singlet and one for the triplet initial spin state, multiply each of 
them by a different enhancement factor. The method we used to do so is described in Appendix 

Would one need to include higher partial waves, one has to compute Cij^i'ji (0) coefficients for 



odd partial waves, add a centrifugal term to the equations ( |3.6D and (3.9) and modify the expression 
for the enhancement, see p^ . 

^This form is most suitable for numerical solutions, while often it is presented not in terms of the CM momentum, 
but rather relative velocity of incoming particles, being equal to v = plm!^ . 
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4. Relic density computation 



The multi-state scenario we have just introduced for the Sommerfeld enhancement is a typical setup 
in which the computation of the thermal relic density for the lightest of such states is strongly 
affected by co-annihilations To model this effect, we will follow here the approach of 

Ref. appropriately modified to include the Sommerfeld enhancements; such treatment will not 
involve approximations in the way thermally averaged annihilation cross sections are computed, 
while the density evolution equation is solved fully numerically. We summarize briefly here the 
steps involved in the computation. 

Let xiy X2, ••• Xm, be a set of M particles, each with mass (the ordering is such that 
mi < 7712 < ■ • • < ™m) and number of internal degrees of freedom /i^, sharing a conserved quantum 
number so that: i) if kinematically allowed, inelastic scatterings on SM thermal bath particles 
can turn each of these states into another, and ii) xi is stable. If the mass splitting between the 
heaviest and the lightest is comparable with mi/20, roughly speaking the freeze out temperature for 
a WIMP, all these states have comparable number densities at decoupling and actively participate 
to the thermal freeze out. A system of M coupled Boltzmann equations can be written to trace 
the number density Ui of each single state; since, however, after freeze out all heavier states decay 
into the lightest, one usually solves a single equation written for the sum of the number densities, 
n = Y^i^r^ i-e- 11: 

5-K3i77i = ~(aeff«) , (4.1) 

at 

with the effective thermally averaged annihilation cross section: 

{a,^y)^J2{a,,v,^)-±.-L. (4.2) 

written as a weighted sum over the thermally averaged annihilation cross section for processes of 
the type Xi Xj ^ ^ (in the dilute limit, two-body initial state processes dominate): 

K%>-;^E / ^^^^s\p,+P,-Px)n'iP^)friP,)\A,^x\\ (4.3) 

where the sum here is on the set X of allowed SM final states (normally only two-body final states 
are considered) and d^px /2 Ex stands symbolically for the integration over the phase space in the 
final state. In the equations above = J (Ppf^'^ (pi) is the thermal equilibrium number density 
for the species i, while n'^'^ the sum of n^'' over all states. There are two main assumptions which 
allow to rewrite the system of M coupled Boltzmann equations as the single equation for n, in 
analogy to the case one writes for one single WIMP, with the usual term of dilution by volume on 
the l.h.s. and the depletion and replenish terms on the r.h.s. The factorization of the individual 



terms in the sum of Eq. (4.2) is possible if one assumes that the shape of phase space densities for 
each particle Xi follows the shape of the corresponding thermal equilibrium phase space density, 
namely /i(pi, t) — c(t) ■ f^''{pi, t) (with the coefficient c depending on time but not on momentum); 
this is the case if the so-called kinetic equilibrium is maintained, i.e. if scattering processes of 
the kind Xi + ^ Xk + Xm (with k equal to i or different from it) on SM thermal bath states 
Xi and Xm have a rate which is larger than the Universe expansion rate H . The amplitudes 
for scattering and annihilation processes are usually comparable since the two can be related via 
crossing symmetry; on the other hand scatterings are stimulated by thermal bath states, which are 
relativistic and hence whose number density is much larger than the number density of the particle 
X at the freeze-out temperature ~ mi/20. Hence kinetic decoupling usually takes place much later 
than chemical freeze-out (chemistry here indicates number density changing processes). Since, as 
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we will show below, within the particle physics setup we are considering the Sommerfeld factors are 
of order 0(1), this is typically the case even when the depletion term takes the Sommerfeld effect 
into account. In case of a resonance the kinetic decoupling can indeed happen when the depletion 
term is still active, as we will discuss in Section ^.5| . The kinetic decoupling can have a much larger 
impact in more extreme scenarios, as discussed in Ref. [Q. 

The second assumption which is implicit in Eq. (4.1) is that one takes n.i/n ~ n^'/n'^'^, a 



quantity which in the Maxwell-Boltzmann approximation for equilibrium phase space densities, as 
appropriate for WIMPs, i.e. fi'^{pi,t) = hi/ {2Tr)^ exp{—Ei/T), is proportional to the number of 
internal degrees of freedom hi and is exponentially suppressed with the mass splittings. Analogously 
to the first assumption, this approximation is valid in case inelastic scatterings are active for the 
whole phase in which the depletion term is relevant. 



Usually the thermally averaged annihilation cross sections in Eq. (4.3) are computed at the 
lowest order in perturbation theory taking tree-level amplitudes. Here we will include the Sommer- 
feld eff'ect introducing the rescaling \Aab\^ — J2ij '^{Lb) l^ijl^' '^ith '^[ib) ^ computed in previous 
Section. Actually, since the effect can be interpreted as a rescaling in the wave function of the 
incoming pair, •S'^^^-j does not depend on the final state X and can be factorizes out of the total 
annihilation rate Wij. Following the same steps of Ref. and adopting an analogous notation, 
one finds: 



[CTcSV) 



(4.4) 



where we have defined: 



5VF,ff(peff,r) = ^^!gi^(^^^^(p,ff,r)w^„(peff) (4.5) 

ab ^ 



with: 



1 / [s - (m, - mjf] [s - (m, + m^Y] 
Pr3 = 2 V 1 (4.6) 

and PcS = Pii = s — Am\. The explicit dependence of S^^^^-^ on T stems from the fact that 

there may be an explicit dependence of the mass of the long-range force carriers on temperature, 
see the discussion below; its dependence on the relative velocity of the annihilating pairs has been 
rewritten instead in term of the integration variable, i.e. the effective momentum poff ■ 



5. Implementation of the method to the MSSM 

We apply the general method outlined in the previous Sections to one of the best motivated and 
most extensively studied scenarios, namely the case of neutralino dark matter in the MSSM. Since 
we will focus the discuss to cases in which the sfermion sector is not playing any relevant role, to 
simplify the discussion and underline better which are the key parameters, we choose actually to 
refer to the SUSY framework usually dubbed "Split Supersymmetry" This indicates a 

generic realization of the SUSY extension to the SM where fermionic superpartners feature a low 
mass spectrum (say at the TeV scale or lower), while scalar superpartners are heavy, with a mass 



scale which can in principle range from hundreds of TeV up to the GUT or the Planck scale [33 
a feature which can occur in wide class of theories, see, e.g. 36, We will also leave out of 
our discussion the Gluino and the gravitino, supposing they are (moderately) heavy, focussing the 
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analysis on neutralinos and charginos.'^ The MSSM contains four neutralinos, spin 1/2 Majorana 
fermions arising as the mass eigenstates from the superposition of the two neutral Gauginos, the 
Bino B and the Wino W^, and the two neutral Higgsino fields Hi and ff^ ■ 

X° = N,iB + 7V,2M^' + N^sH^i + N,^H^ ; (5.1) 

in the scheme we are considering the lightest neutralino is automatically the lightest SUSY particle 
(LSP) and, possibly, a good WIMP dark matter candidate. The coefficients Nij, obtained by 
diagonalizing the neutralino mass matrix, are mainly a function of the Bino and the Wino mass 
parameters A'li and M2, and of the Higgs superfield parameter /i, while depend rather weakly on 
tan /?, the ratio of the vacuum expectation values of the two neutral components of the SU(2) 
Higgs doublets, which appears in the off-diagonal terms of neutralino mass matrix. The two MSSM 
charginos are instead spin 1/2 Dirac fermions obtained as mass eigenstates from one charged Wino 
state and one charged Higgsino; the chargino composition is again mainly set by the relative weight 
of M2 and ^. 

The long-range forces driving a sizable Sommerfeld effect are the weak force on neutralinos and 
charginos due to the exchange of and 2'°, and the electromagnetic force with photon exchanges; 
we also include the effect of the charged Higgs and the two CP-even neutral Higgs Hi and H2, 
which however play a minor role. The last MSSM scalar, the CP-odd neutral Higgs H^, does not 
give rise to any contribution for s-wave annihilations in the non-relativistic regime. Since the B is 
not charged under SU{2)l, a large Bino component in the lightest neutralino drastically reduces the 
relevance of the Sommerfeld effect; we will then consider only the case of Mi ^ M2. Higgsinos and 
Winos have a pair annihilation cross section into W and Z bosons which is fairly large, much larger 
than the standard reference value for thermal relics of 3 • 10^^^ cm"^ s^^ if their mass is at around 
100 GeV. Going however for more massive neutralinos, namely around 1.1 TeV for Higgsinos and 
2.2 TeV for Winos, the standard tree-level calculation of the thermal relic density gives a result 
which is compatible with the measured value for the energy density in CDM. This heavy mass 
regime is also the one in which the Sommerfeld enhancement condition of mass of the particle much 
heavier than mass of the force carrier is realized for weak interactions; hence relevant corrections to 
the tree-level estimate of the relic abundance may arise |l|,|l|. Ah the pair annihilation processes 
we will need to consider are dominated by their s-wave contribution (unlike the case of the pair 
annihilation of Binos into a fermion-antifermion) ; moreover the enhancement for the higher partial 
waves is smaller, due to repulsive centrifugal term in the Schrodinger equation.^ We can then safely 
assume that only the s-wave Sommerfeld effect is relevant. 

As we stressed in Section ^, the value of the mass splitting among the different states is crucial 
in the analysis. When calculating the mass spectrum we include the radiative corrections to the 
neutralino and chargino masses due to gauge boson loops |^ , 40 (tree-level neutralino and chargino 



masses are degenerate in the pure Higgsino or pure Wino limits) and the thermal dependence of the 
Higgs VEVs. Thermal effects are also important when considering vector boson masses, as already 
stressed in the analysis of Ref. [O]. There are two thermal effects which we need to include: First, 



''The Sommerfeld effect does play a role in the pair annihilation rate of charged sfermions or gluinos, and so it 
does affect the computation of the neutralino relic density in case of sfermion or gluino co-annihilations; there is 
however no long-range interaction turning a pair of these particle into a neutralino pair, and hence, contrary to case 
discussed in this paper, the different annihilation can be treated separately and the calculation is simplified. An 
analysis dedicated to these co-annihilation channels will be presented in Ref. 

■'This statement is not true only in the case of resonances. In the Coulomb case, when analytical expressions can 
be derived, every partial wave is enhanced by the factor ~ (a/v)^'^^ , so for higher I the resonant enhancement is 
higher and narrower. However, since one has to integrate over the thermal distributions, the resonances are smeared 
and net impact of the higher partial waves is very small. 
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we approximate the scaling of the VEV v with temperature as 



viT) = vRe (^1 - ^2 j , (5.2) 

with the critical temperature Tc depending on the Higgs mass (as in Ref. [|2j, we will assume 
Tc = 200 GeV). Second, we consider the contribution to gauge boson masses, or, in more appropriate 
terms, to their propagator pole, due to the screening by the thermal plasma; this effect can be 
approximated by adding the so-called Debye mass H), which in the SM and at T 'rnw,z are: 

We will assume that these expressions are valid also for the MSSM and at lower temperatures; the 
first assumption is justified by the fact that the contributions from the additional states are small 
since these particles are heavy, the latter has a negligible impact on the relic density calculation. 
The two effects introduce a correction to the Sommerfeld factors which is non- negligible, but still 
quite small, since they are important mostly for high temperatures, at which the Sommerfeld effect 
is anyway negligible. 

Having computed with high accuracy the mass spectra, we implement the procedure for the relic 
density, considering a system of coupled equations which includes all states with small mass splitting 
compared to the LSP. In practice this reduced to including two Majorana and one Dirac fermion 
when fi < M2 (i.e. two neutralinos and one chargino, mainly Higgsino-like) and one Majorana 
and one Dirac fermion when /i > M2 (i.e. one neutralino and one chargino, mainly Wino-like). 
Additionally, since the Sommerfeld enhancements depend on the total spin of the initial pair, we 
considered each case separately if needed. Again, in practice this is important only if the incoming 
particles are one neutralino and one chargino: two identical Majoranas cannot form a s-wave spin 
triplet state, the triplet of two different Majoranas has suppressed annihilations (it cannot annihilate 
to M^^W^~), and for two charginos the effect comes dominantly via 7 exchange, which is vector and 
both singlet and triplet computations coincide. 

6. Results and discussion 

Having introduced the method and particle physics framework, we present results for the Sommer- 
feld factors and the impact of the effect on the relic density of the neutralino. As already explained, 
M2 and /i are the critical parameter in this problem; we will consider an approach with parameters 
set at the low energy scale (rather than at a GUT scale as often done) and let them vary freely. 
The other MSSM parameters are kept fixed: besides the sfermion sector which is assumed to be 
heavy, and the Bino which is also decoupled with the artifact of setting Mi = IOOM2, we assume 
tan /3 = 30 and an Higgs sector with a light SM-like Higgs and all other states that are heavy, as 
expected in split SUSY. The value of tan /3 has a very modest impact on results and considering 
other values does not bring other information. In the same way, the other Higgs sector parameters 
have very little relevance. For any given point in this two-dimensional parameter space, the relic 
density calculation is performed computing first the coefficients Sij{v,T) numerically, adopting a 
standard adaptive 5**^ order Runge-Kutta method with shooting for solving the boundary value 
problem [||]. 

6.1 Sommerfeld factors and the effective cross section 

For any given point in the MSSM parameter space leading to a lightest neutralino x° which is 
heavy and has a large Wino or Higgsino fraction there are several Sommerfeld enhancement factors 
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^}ab) useded for a relic density computation. The interplay between the different contributions 

to the thermal averaged effective cross section is in most cases non-trivial. In Fig. ^ we give an 
example of '5'(^f,)(w, T) for the case when the neutralino is nearly purely W^, just to give a intuition 
of what magnitude of enhancement factors we deal with. In this case in set of coupled equations 
(3.9) we have N < 2, i.e. only neutralino and chargino coupled together. One can see that the 
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Figure 2: Example enhancement factors for M2 — 2.53 TeV, fi = 5 TeV and x — 20 (corresponding to 
T ~ 126 GeV). Subscripts refer to the incoming state with 1 - Xi s-nd ± - x^i while superscripts to the 
annihilating one. The left panel shows Sommerfeld factors when the annihilating pair is the same as the 
incoming one, while right panel when they are different. In the latter case these factors approach because 
at high momenta there is no Sommerfeld enhancement and the annihilation amplitude is suppressed because 
it can only be obtained at 1-loop level at least. All factors are computed for singlet spin state annihilation, 
except one indicated as tS being for triplet. 



enhancements are of order 0(1) and very quickly go to 1 or to in higher momenta and that most 
channels are attractive while two are repulsive, namely x~^Xi in ths singlet spin state and x^X^- 
Also, as expected, one gets the highest enhancement factor for the x^X~ channel, since there is 
a long range nearly^ Coulomb interaction present. For the cases with two coupled equations and 
lighter states incoming we see a resonance in the value of momentum corresponding to the mass 
splitting. This is easy to understand, since for this energy the heavier states are produced nearly 
at rest, so they feel large effect coming from for instance 7 exchange. This is perfectly consistent 
with what was found in similar cases using analytical approximations . 

We are now ready to compute Ccff and solve the Boltzmann equation as discussed in Section ^. 
Values for the effective annihilation cross section including the Sommerfeld effect and without it 
are shown in Fig. |^ for sample cases of Wino- and Higgsino-like neutralino. 

In the Higgsino-like case the effect is very mild. The Sommerfeld enhancement of the effective 
cross section becomes relevant only in the small velocities regime, when the depletion term in the 
Boltzmann equation is marginally effective. This gives rise to a change in the relic density which is 
at most at the level of a few per cent. In the Wino-like case the picture looks much more interesting. 
The net effect on the both yield and (Toff is clearly visible and can become even very large in the 
parameter range where large resonance effects occur. 

The reason for the difference in the behaviour of the Sommerfeld effect in the Wino- and 
Higgsino-like case comes mainly from the mass splitting between the lightest neutralino and chargino, 
which is typically much smaller in the Wino case. The Sommerfeld effect for the neutralinos re- 
lies mostly on a production of nearly on-shell charginos in the loop. Also the efficiency of co- 
annihilation, and subsequently the effect of the Sommerfeld enhancement coming from its impact 



^Due to the thermal corrections photon acquires small mass, which makes the potential to be Yukawa type, and 
the enhancement factor saturates at small p. 
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Figure 3: Number densities in units of entropy density Y and effective cross sections for tiie Wino-like (to-p) 
and Higgsino-like neutralino (bottom). In the Wino case two set of parameters are presented - generic one 
and close to the resonance. 



on co-annihilating particles, strongly depends on the mass splitting. Hence, the larger the mass 
splitting the smaller the overall impact on the thermally averaged cross section.^ 

6.2 Relic density 

In Fig. ^ and ^ we show results for the neutralino relic density in the plane M2 versus fx, varying 
these parameters in the range 500 GeV to 5 TeV, a region in which the thermal relic density 
varies from values much below the cosmologically preferred one to much above it. In the left 
panel of Fig. ^, we are plotting results assuming the usual tree-level approximation for the pair 
annihilation amplitude, while in the right panel those including the full treatment of the Sommerfeld 
effect are shown. Most manifestly, there is a sharp shift in the Wino-like region consistent with 
the 7-year WMAP data to heavier masses; when including the Sommerfeld enhancement a pure 
Wino is found to have fth^ = 0.11 for a mass of about 2.8 TeV. Much milder changes place in 
the Higgsino-like region; the relic density is practically unchanged when considering models in the 
cosmologically interesting band. The relic density decrease is even larger in the ultra-heavy regime, 
with the Sommerfeld effect becoming larger and larger as the gauge boson masses are becoming 
less important; this regime would however be consistent with cosmology only invoking some extra 
ingredient, such as, e.g., a dilution effect via late entropy production or the decay of the neutralino 
into a lighter state, such as a gravitino or an axino, which is the true LSP (although viable, both 
scenario require large fine-tuning). The results found here for pure Winos and pure Higgsino are 
analogous to those of earlier works in Ref. and Ref. [|l2| in the same limiting cases; there are 
small quantitative differences stemming from the fact that we have identified more annihilation 
channels that need to be treated separately, with the Sommerfeld factor depending on the initial 
spin state of the annihilating particle pair, we have found a different coefficient in the axial vector 

^Note however, that in general Sommerfeld factors themselves are not monotonic functions of the mass spUtting, 
see e.g. Ref. Q. 
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exchange, i.e. that the axial vector has an additional —3 factor with respect to the vector (in 
agreement with the result in Ref. and we have probably a better control of the numerical 
solution of the Boltzmann Eq. having implemented our full treatment in the DarkSUSY numerical 
package (the slight difference in numerical results between [pi and are instead probably mainly 
due to the thermal corrections implemented here following [f2|). 
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Figure 4: Relic density 0,h^ in the fi-M2 plane for perturbative case (left panel) and with Sommerfeld effect 
included (right panel). The brighter the colour the higher Q,h^ and the colour scale is linear. The solid line 
and dashed lines correspond to the central value and the 1 a error bar for relic density consistent with the 
7-year WMAP data, Qh'^ = 0.1123 ± 0.0035. 




Another region showing interesting results is the band connecting the pure Higgsino to the pure 
Wino limit, towards M2 ~ /i but still with a predominant Wino component. Features in this region 
are more clearly seen in Fig. ^ where we show the ratio between the relic density computed with 
tree-level amplitudes to the one with the full non-perturbative treatment. A thin "resonance" slice 
appears in the plane, starting for pure Winos with mass m^o « 2.5 TeV and extending to heavier 
masses into the region with a sizable Higgsino fraction, where the thermal relic density becomes 
consistent with observations. The value of the mass we find for a pure Wino is precisely the one 
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saturating Eq. (2.5), i.e 



(6.1) 



with a as computed in the vertex for pure Winos W'^ W~ . This means that the observed 
resonance is due to the possibility of creating the loosely bound state, occurring when the Bohr 
radius coincides with the interaction range. It also explains why, when we increase the Higgsino 
fraction, the resonance deviates to higher masses: a larger Higgsino fraction implies an increase in 
the mass splitting between lightest neutralino and chargino (since one goes from a mass splitting 
dominated by the radiative corrections, about 170 MeV, to the one induced by the mixing of 
interaction eigenstates) , as well as a drop in the couplings (since the vertex iJ" i7+ W~ as a coupling 
which is a factor of \/2 smaller than for Winos) and hence a drop in a; this has to be compensated 
by a larger m^o . 

6.3 Kinetic decoupling 

The results we presented in previous Subsections assume that the temperature of the neutralinos 
traces the thermal bath temperature. This is true as long as neutralinos are in kinetic equilibrium. 
After kinetic decoupling, at the temperature Tkd, their temperature decreases with the scaling as 
appropriate for non-relativistic particles, i.e. ~ 1/a^, where a is the Universe scale factor, 
cooling much faster than thermal bath states for which T ^ 1/a. This does not have any influence 
on the relic density computation if interactions do not depend on velocity, as in the standard 
case of s-wave annihilations. When the Sommerfeld effect plays a major role, however, there is 
indeed a strong dependence on velocity and colder neutralino may have larger annihilation cross 
section. Hence, if the kinetic decoupling happened early enough, i.e. when the depletion term in the 
Boltzmann equation is still active, this might have given rise to stronger relic density suppression 
by the Sommerfeld enhancement. 
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Figure 6: Influence of the kinetic decoupling on the relic density suppression by the Sommerfeld efi'ect for 
Wino-like neutralinos. In all cases fj, = 5 TeV and other parameters are as before. For a pure W case the 
kinetic decoupling occurs at about Tkd ~ 170 MeV, as explained in the text. There is a clear enhancement 
in the net effect within the resonance region, but no effect away from it. 



We have just shown that, within the MSSM, the largest impact of the Sommerfeld effect occurs 
for Wino-like neutralino. In this case, since elastic scatterings are very strongly suppressed, the 
processes enforcing kinetic equilibrium are inelastic scatterings of the type Xi^~ ^ X~Ve\ the latter 
are efRcient up to the time when the temperature drops down below the mass splitting between 
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neutralino and chargino. Hence, we can estimate that the kinetic decouphng for Winos occurs 
at about Tkd ~ 170 MeV. The exact value of Tkd and shape of the distribution function after 
decouphng (since it's unhkely for the decouphng to be instantaneous) depend on the parameters 
in the model; it is in general rather involved to determine them (see e.g. ^) and this goes 
beyond the purpose of this paper. Here, to illustrate the possible effect of the kinetic decouphng on 
the relic density when the Sommerfeld enhancement is relevant, in Fig. ^we plot the ratio between 
the value of the relic density as computed at the lowest order in perturbation theory (i7/i^)o and 
the full computation (r2/i^)gE, for a few points in the parameter space with Wino-like neutralinos, 
as a function of the value assumed for Tkd and assuming the neutralino distribution function as 
if decoupling is instantaneous. One can see that in general, the relic density is not sensitive to 
the kinetic decoupling temperature; the exception is the resonance case, when the Sommerfeld 
enhancement can be much larger and we find a sizable corrections depending on Tkd. It follows 
that accurate predictions of the relic abundance in the resonance regime are possible only after 
determining the kinetic decoupling temperature with a certain accuracy. On the other hand the 
resonance region is rather tiny and the overall MSSM picture discussed in this work is not much 
affected. 

7. Conclusions 

In this paper we have computed the Sommerfeld effect in neutralino and chargino pair annihilations 
within the MSSM, and discussed its relevance for the thermal relic density of the lightest neutralino. 
For this purpose, we have adopted and extended the general approach presented in to obtain 
a reliable method of computing the Sommerfeld factors in a case of multi-state systems of fermions 
with interactions mediated by scalar, vector and axial vector bosons. 

The same method is readily extendable to other models independently of the number of states 
and interactions that are present; e.g. one could use this same approach within next-to minimal 
SUSY models (NMSSM), where in some cases a large Sommerfeld effect may be expected due to 
presence of an additional chiral superfield coupled to Higgs fields. Moreover, the results presented in 
this analysis are already readily readable for more general models than the MSSM: the ingredients 
we have implemented are only: i) the existence of massive fermions in the adjoint representation 
of SU{2)l (we have been calling them Winos); ii) the presence of massive left-handed fermions 
in the fundamental of SU(2)l (without further right-handed ones), which form two doublets with 
hypercharge assignments such that a SU (2) ^ gauge invariant mass is possible (like the /i-term; we 
have been calling these Higgsinos); in) they form the dark matter and that their SU{2)l gauge 
invariant masses are in the range of TeV; iv) there are SU{2)l breaking corrections to the mass 
matrix in the off-diagonal terms in the range of the EW breaking scale, or also in the diagonal 
terms with a scale comparable to the effect of the off-diagonal ones, such that the resulting masses 
receive a correction of the order of O(10~^) with respect to their SU{2)l invariant masses (these 
determine the mixing among neutral and among charged states; we have been referring to the 
mass eigenstates from the superposition of interaction eigenstates as, respectively, neutralinos and 
charginos). Regarding the last ingredient, we have been computing mass matrices and radiative 
corrections under the standard MSSM assumptions but this should be just taken as a specific 
parameter fixing in a more generic model which can be readily generalized. 

For our sample numerical study, we have concentrated on neutralino dark matter in the MSSM, 
being one of the best motivated and most extensively studied scenarios. Since the Sommerfeld effect 
is relevant only for non-trivial representations of SU{2)l we have assumed that the Bino is much 
heavier and considered only the case of Wino-Higgsino mixing; also since in this limit sfermions play 
a marginal role, we have restricted the analysis to a split SUSY scenario, making the hypothesis 
that all scalars except for the SM-like Higgs are heavy. The computation of the relic density is 
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based on a numerical solution of the system of coupled Schrodinger equations which allows to 
compute Sommerfeld factors in all (co-)annihilation channels and then on their implementation in 
the DarkSUSY package for a full numerical solution of the Boltzmann equation. This enables us to 
determine the thermal relic density with a very good accuracy, possibly at the few per cent level, in 
all the parameter space, except for a small region in which a resonance is found because the Bohr 
radius becomes of the same order as the interaction range and for which a careful computation of 
the kinetic decoupling process would be needed. 

The results we have found show that the Sommerfeld effect can suppress the relic density of 
the neutralino by as much as a factor of 2-3 in sizable regions of parameter space and even much 
more in the resonance (the largest suppression factor found in our scan was about 6.2). If the 
LSP is mostly Wino this results in the shift of the neutralino mass being compatible with WMAP 
from about 2.2 TeV to about 2.8 TeV, also enlarging slightly the region allowed under standard 
assumptions on the cosmological model. The existence of the resonance has a small, but visible 
impact on this regime preferred by cosmology, however it may be much more important if one does 
not require the neutralino to constitute the whole CDM term. Similar features may occur in other 
regions of the parameter space; we defer the study of these to future work. 

Qualitatively, the results we have found for pure Winos and Higgsinos are analogous to those 
in earlier studies of these models, see |l^; quantitatively, they slightly differ. There are several 
reasons for that. First of all we were computing the relic density with a full numerical treatment, 
properly including co-annihilations and Sommerfeld factors, having implemented our new treatment 
into the DarkSUSY package. Secondly, we have revised the Sommerfeld enhancement computation 
itself: we have identified more annihilation channels that need to be treated separately since we 
have noticed that the enhancement factor depends on the initial spin state of annihilating particle 
pair, and we have corrected one of the couplings. We also included Higgses exchanges and treated 
separately two different neutralino states (however, this latter effect is quantitatively not relevant 
within our parameter space). Moreover, contrary to previous results, our approach allows us to 
treat not only pure states, but the more realistic situation with arbitrary Higgsino-Wino mixing, 
with mass spectrum as computed in an arbitrary setup, hence allowing us to explore the full MSSM 
parameter space. 

As a final remark, we would like to point out that, although the numerical work carried out 
in this paper has been concentrated on the impact on the relic density calculation, the method we 
developed and the relative numerical routines can be also used to give more accurate predictions 
for the pair annihilation rate of neutralinos in dark matter halos today and hence for indirect dark 
matter detection signals; we will investigate this further in future work. 
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A. Method of computation 

The form of the (r) for any two-particle states exchanging a boson cj) is always of Yukawa- or 

Coulomb-type, but due to different couplings and multiplicities we can have different relative coef- 
ficients in front. Here we will summarize the way of computing them. Those numerical coefficients 
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depend on the fact what type of fermions, Dirac or Majorana, are present in the diagram and are 
they distinguishable or not. 

In it was shown that the integration on the time component of the loop momentum of the 
Feynman graph expression of the Bethe-Salpeter equation for a four-point amplitude gives in 
non-relativistic limit: 

Ha - c 

From this expression, by redefining ^ = {Hq — we get the Schrodinger equation 

i?o$ - V,nt = (A.2) 

Here $ is in general a multicomponent state describing many possible pairs of particles that interact 
with each other, therefore the equation is just a compact form for a system of equations, and the 
symbol * states for an operator acting on spins as well. For each pair ij: 



2mr 



^o' =-:r^ + 2'5m,,, (A.3) 



where my = rriimj / {mi + rrij) is the reduced mass of the pair and 2Smij = mi + mj — {nia + mfj) 
where ma and m^ are the masses of the incoming particles. £ is the kinetic energy (at infinity) of 
the incoming pair. 

In order to make explicit the action of the interaction Vint on a state, we express both the 
interaction and the state in terms of fields and then make the appropriate contractions. The 
relativistic Feynman diagram in the non-relativistic limit gives explicitly the result, but it is simpler 
to work directly with the form ( [A.2| ) by defining suitable non-relativistic contractions. This way 
of doing the computation has a virtue of easy bookkeeping of all multiplicative factors and signs 
appearing especially in the case involving Majorana particles. This formalism allows also to include 
easily the fact that we have initial and final state with defined total angular momentum and spin. 

In the non-relativistic approximation the time delay is neglected and we can work in the time 
independent Schrodinger picture. 

We introduce the state $ describing a fermion-antifermion (Dirac or Majorana) pair expressed 

as' 



1$;^) - j dzdwMz)O^^jiw)\0)%'iz,w). (A.4) 

For a (Dirac) fermion-fermion pair one needs to take iplw) — > ip'^{w). It is easy to see that the spin 
singlet 5 = and the spin triplet S* = 1 are encoded in this formula*: 

S^Q: 0^ = 75, 5=1; 0^ = ^-S, (A.5) 

where S is the spin of initial pair. The normalization is 

N,, = I/V2 z ^ J , (A.6) 
N^, = 1/2 i = j. 

Take an interaction of the form (as usual by F we denote the gamma matrices structure) 

Vnt =9l I dxdyMx)Ti^^{x)i^M^Mv)wtu,ix - V). (A.7) 



^In the following we write x, y, z, w for x, p . z, w, the time coordinate (not indicated) being the same. 
*Here we extend the idea presented in ]45| to include also the spin triplet. 
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where W^j^i is the propagator of the boson exchanged between the two vertices. In the non- 
relativistic hmit only T = 1,70,7^75 can contribute. The transition ij — >■ kl can be described in 
terms of operators acting on the initial state giving the final one: 

|i J dxdy {M^)T^|^^{x) + Ux)Ti;,{x) + h.c.) {My)^My) + My)^i'iiy) + h.c) W^^^^ilx ~ y\)x 
X j dzdv3Niji!,{z)0^i;j{w)\Q)^'^[z,w) ^ 

dxdyNMMx)O^MymV^i^^^{\x-y\)^'^{x,y). (A.8) 

The interaction potential between two 2-particle states {ij — >■ kl) arising due to vector, axial vector 
boson or scalar exchange with the mass has the form (see |26| for details on the normalization) 

y.Ur)^^'—-. (A.9) 

We are interested in computing the coefficients cj^^ ^^^^-^ for all possible cases. Note that in the non- 
relativistic approximation spin-orbit interactions are suppressed, therefore a spin-singiet (triplet) 
initial state gives a a spin-singlet(triplet) final state, and also the parity of the wave function, 
that is whether ^(r) — ±<I>(— r) (e.g. + for the s-wave and — for the p-wave), is the same in the 
initial and final state. Therefore in general there are four independent systems of equations: the 
spin-singlet even, spin-singlet odd, spin-triplet even, spin-triplet odd. The coefficients c^j and 
therefore the interaction potentials are in general different in the four cases. 

To illustrate the method, consider first a (Dirac) fermion-antifermion pair. We make the possible 
contractions: 

v,nt * 1$!/) = 5' / Mx)r{i^,{x)Mz))o^{i'jM^jiy))^My)wHx - y)%'{z,w). 

J xyzw 

By noting that only the creation operator part of both ijja{z) and ^b{w) appear in the state 
il)i{z)0^'ipj{'w)\Q), we get for p <^m: 

= {Omx)M^m = J ^-^e^?(--T ^ usu, ^ S{x-z)P+ , (A.IO) 

s 



where P± = and w = ^ (p- + m^. Therefore 

y^nt * l^-;^) = .g' / dxdyi;kix)rp+0-,p^rMymw^i,^j{^ - yW^{x,y). 



Note that Mx)P+\0) = M^M and that P^MvM = MyM- 

Keeping into account the sign difference between the vector and axial propagator with respect 
to scalar one and defining W{r) = e~"^'t'^ / 4:7rr , one can write W-y^ = —W, Wi^-y^-y^ = +W. Then by 
Dirac algebra we finally get 

V,nt * 1$!/) = y dxdyM=c)O^MyMWix~y)^\\x,y), (A.ll) 

with 

„1 _ „7o _ „1 _ „7o _ 1 ^7375 _ _o ^7j75 _ 1 

75 7i 7i ' 75 ' 7i 
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This result gives a term in the equation for 



where c^j ^^- = j^<^^g^. 

The state may contain Majorana fermions x — x'^ = C'x^j like 

i,{z)0^x{w)\m{^. w). x{z)0^^{w)\Qmz, w), x»(2)O^Xj (w)|0)$(^, w). 

Again = 75 for the spin singlet and 0~, = 7j for the spin triplet. Since x{^) ^-iid x(w) contain 
the creation operator part only, the possible contractions between Majorana fermions are, in the 
non-relativistic limit, 

{x{x)x{z)) ^ 5{x ~ z)P+ , {x{w)x{v)) 5{y - w)P^ , (A.13) 
{X{x)x{wf) = {x{^)xH)C^ ^ ~S{x-w)P+C, 

ixiyfxi^)) = C^{x{v)m) ^ -5{y- z)CP+ . 

The computations presented above are for the even orbital angular momentum, for instance I = 
(s-wave). This means that the 2-body wave function is symmetric. The results for these coefRcients 
are summarized in Table |^. It is very easy to generalize it to higher partial waves. 

These results cover all the possibilities of types of 2-body to 2-body interactions when we can 
have Dirac or Majorana states, i.e. every other case present in some model will have one of those 
forms (so the coefficient in Schrodinger equation for computing Sommerfeld enhancement will be 
the same) with possible different coupling constant. 

B. Projection of the (co-) annihilation amplitudes into spin singlet and 
spin triplet initial states 

Since the Sommerfeld effect is in general different for the spin singlet and triplet configurations, it 
will affect separately the annihilation rates for those two cases. Therefore, we need to compute the 
relative weight of the two cases. 

We consider the annihilation of two fermionic dark matter particles into two bosonic or fermionic 
standard model particles. In our approximation we take the initial state to be at rest, and we neglect 
the dark matter mass differences and the masses of the particles resulting from the annihilation. 
Hence the kinematics of the annihilation process simplifies; by calling Pfi,p'^ the momenta of the 
incoming dark matter particles and fc^ , k'^ the momenta of the annihilation products we get in the 
CM frame^: 

Pp--(m, 0), p'^^{m,Q), k^j_ ^ {'m,mn) , k'^^{m,-mn), = I . (B.l) 

The propagators of the virtual particles exchanged in the Feynman diagrams representing the 
amplitude will have denominators of the kind: 

t-m^ ^ -2m^ , u- ^ -2m^ , s - fi^ ^ 4m^ , (B.2) 

where t = {p— k)"^ ~ —m?, u= {p — k')^ ~ —m^, s={p + p'Y ^ 4m^, and m is the appropriate 
mass of the particle mediating the annihilation. 

^It can be done without this simphfication but in the case of interest it does not play any role, especially that 
we are interested only in the relative weights and so by doing this approximation we do not change the value of 
perturbative annihilation cross section. 
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With these simplifications one can decompose the rate to two factors, i.e. rate (x K x ^, where 
K is the amplitude squared summed over the final spin configurations and averaged over the initial 
one, and $ is the phase space. 

The phase space factor, apart from common factors, reads 



C7n 



cm 



E^ 

cm 



1/2 



(B.3) 



where mi and m2 are masses of final particles and £'^„j is the energy in the CM frame 




Figure 7: Diagram with explicit contraction of the initial spinors by O-y. The blob represents every possible 
annihilation process with all possible 2-body final states - fermionic or bosonic. 



The relative weights for singlet (triplet) will be computed by the ratio 



Q(5 = 0,l) 



Q{S = 0) + Q{S ^ 1) ' 

where Q{S = 0, 1) is the sum of the singlet (triplet) squared amplitudes for the various annihilation 
channels: 

Q{S = 0, 1) = ^ Qchannel j{S = 0, 1). 
j 

In conclusion, in our approximation, the Sommerfeld enhanced total (i.e. summed over every 
channel) annihilation cross section aenh is related to the non-enhanced total annihilation cross 
section (Tq by the formula 



^enh 



Q{S = 0) + Q{S = I) 



S{S = 0) 



Q{S = 1) 



Q(5 = 0) + Q(5 = 1 



■5(5 = 1) 



(B.4) 



where S{S = 0, 1) are the Sommerfeld enhancement factors, proportional to the square modulus of 
the non-relativistic wave function at the origin, for the spin singlet and triplet. 

The computation of the weights follows closely the method previously seen. We compute the 
standard tree-level amplitude of annihilation to given final state but with initial spinors contracted 
by O-),, as presented with a diagram on a Fig. ??. In a non-relativistic limit this gives precisely 
amplitudes Q{S = 0) for 75 and Q{S — 1) for 7 • S. From that point the computations proceed in 
a standard way. 
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